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1. Introduction 



Asymptotically free theories, such as QCD, behave very differently at low 
and high energies. While there is no obvious physical relation between, say, 
pion-nucleon scattering phases and jet production rates at e"'"e~ colliders, the 
underlying theory describes both phenomena and thus provides a link between 
the two. 

An interesting observation in this context is that the running coupling at 
high energies may in principle be computed once the parameters of the theory 
are fixed at low energies [1]. This computation is difficult, because one refers 
to the non-perturbative low-energy properties of the theory. Moreover, one 
must be able to cover a substantial range of energies to make contact with 
the scaling regime, where the running coupling is logarithmically decreasing 
according to the perturbative renormalization group. 

In this paper the problem is solved for the SU(3) Yang-Mills theory. The 
computation is set up in the framework of lattice gauge theories and relies 
on numerical simulations. To step up the energy scale, a recursive finite-size 
scaling technique is employed. The method has previously been applied to the 
SU(2) theory [3,4] and we shall here assume that the reader is familiar with 
this work (for other approaches see refs.[6-ll]). 

There are two new features in our computations. The first is that we 
now use a 1-loop improved lattice action, which leads to significantly reduced 
lattice spacing effects. This makes the extrapolation of the lattice data to the 
continuum limit safer. 

The other change concerns the choice of reference scale at low energies. In 
the SU(2) theory all physical momenta were given in units of the string tension 
K. From a technical point of view this quantity is somewhat problematic, 
because in practice it is determined by extrapolating the heavy quark force 
F[r) from distances r less than 1 fm to large distances. Depending on which 
analytical form is assumed for the extrapolation, one can obtain quite different 
results. In other words, the string tension is affected by a systematic error 
which is not easy to control. 

An alternative reference scale is the distance t-q at which 

rlF{ro) = 1.65 (1.1) 

The merits of this definition have been discussed in ref.[5]. The most important 
points to note are the following. 
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1. ro is well-defined in SU(iV) gauge theories with or without matter fields. 

2. From the non-relativistic charmonium model one estimates that t-q is ap- 
proximately equal to 0.5 fm in nature [the number on the right hand side of 
eq.(l.l) has been chosen so as to achieve this]. The exact relation of t-q to 
other scales in QCD, such as the pion decay constant, is presently not known. 

3. In lattice gauge theories t-q is comparatively easy to compute through 
numerical simulation, because one only requires the static quark force in an 
accessible range of distances. No extrapolation is needed here. 

4. In the SU(3) Yang-Mills theory one finds 



where the error includes an estimate of the systematic uncertainty on the 
string tension mentioned above (cf. sect. 5). 

Neither the distance t-q nor the string tension are experimentally measurable, 
but ro has several conceptual and technical advantages so that we decided to 
express our final results in units of this scale. 



The study reported in this paper is a straightforward extension of our 
earlier work on the SU(2) Yang-Mills theory [2-4]. We assume the reader is 
familiar with our computational strategy. Here we only introduce the neces- 
sary notations. 

We choose to set up the theory on a hyper-cubic euclidean lattice with 
spacing a and size Lx Lx Lx L. The possible values of the time coordinate a;" 
of a lattice point x are a;" = 0, a, 2a, . . .,L [L is taken to be an integer multiple 
of a). The spatial sublattices at fixed times are thought to be wrapped on a 
torus, i.e. we assume periodic boundary conditions in these directions. 

A gauge field U on the lattice is an assignment of a matrix U{x, /x) G SU(3) 
to every pair [x,x + afi) of nearest neighbour lattice points [fi denotes the 
unit vector in the /x-direction and /x = 0,1,2,3). At the top and bottom of 
the lattice, the link variables are required to satisfy inhomogeneous Dirichlet 




(1.2) 



2. Lattice theory 
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boundary conditions, 



U{x, fc)|,o=o = k), U{x, k)lo^^ = W'{k, k), (2.1) 

for all = 1, 2, 3, where W and W' are prescribed spatial gauge fields. They 
will be set to some particular values below. 
The partition function of the system, 

Z= /D[C/]e-^[^l, I)[U] = ]ldU{x,fi), (2.2) 

involves an integration over all fields U with fixed boundary values W and 
W . For the action S[U] we take 

S[U] = \Y.'"{p)tv{l-U{p)}, (2.3) 

% p 

with go being the bare coupling. The sum in eq.(2.3) runs over all oriented 
plaquettes p on the lattice and U{p) denotes the parallel transporter around 
p. The weights w[p) depend on whether one desires to work with the simple 
Wilson action or an improved action. In the first instance w[p) is equal to 1 
except for the spatial plaquettes at a;" = and = L which are given the 
weight i. 

As discussed in sect. 4 of ref.[2], the leading cutoff" effects in the Wilson 
theory are of order a. They can be cancelled by adjusting the weights of the 
plaquettes at the boundary of the lattice, i.e. one sets 

\c.s{go) if p is a spatial plaquette at a;" = or x^ = L, 
Ct{go) if p is a time-like plaquette attached to a 

(2.4) 

boundary plane, 

and w[p) = 1 in all other cases. With appropriately chosen coefficients Cs{go) 
and Ct[go), the cutoff effects are then reduced to order . 

We shall exclusively be concerned with constant abelian boundary fields 
W and W' . Since the contribution of the spatial boundary plaquettes to the 
action vanishes in this case, we do not need to know the value of Cs{go)- For 
the other improvement coefficient we use the 1-loop formula 

ct{go) = 1 + c[^\l = -0.08900(5). (2.5) 



w[p) = < 
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This removes all cutoff effects of order a at the tree and the 1-loop level of 
perturbation theory. 



3. Running coupling 

From the partition function Z (which is also referred to as the Schrodinger 
functional) a running coupling may be defined by differentiating with respect 
to the boundary values [2]. To obtain a unique coupling, the boundary fields 
and the direction of differentiation must be specified. There are no deep 
theoretical reasons for the choices made below. They are, however, based on 
practical considerations, such as the requirement of mild cutoff effects and 
calculability. 

3. 1 Boundary values and induced background field 

As in the SU(2) theory we choose the boundary values W and W to be 
constant diagonal matrices. More precisely we set 

W^(x,fc) = exp{aCfc} and W^'(x, fc) = exp{aC;}, (3.1) 

where 

. /./.I \ . /./-i \ 

Cfc = y h and = I </.^ . (3.2) 

VO ./.a/ ""Xq Q <i>'J 

The angles ^„ and (f)'^ must be real and they should add up to zero to guarantee 
that the boundary link variables are elements of SU(3). 
The classical field 

F(a;,/x) = exp{a5^(a;)}, (3.3) 

with 

5o = 0, Bk=[x^C'k + {L-x^)Ck]lL, (3.4) 

is a solution of the lattice field equations with the required boundary values. 
V is referred to as the induced background field. A theorem proved in ref.[2] 
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asserts that V is absolutely stable if L/a > 5 and if the vectors (j)2, (ps) 
and (f)'2,(t>'^) are in a certain bounded region, the fundamental domain. 

Stability ensures that at small couplings go, the Schrodinger functional 
is dominated by the field configurations in a neighbourhood of the induced 
background field. The perturbation expansion then amounts to a saddle point 
expansion about V. 

3.2 Geometry of the fundamental domain 

The fundamental domain consists of all vectors (^i, V^'s) satisfying 

(t)i + (t)2 + h = 0, (t>i <(t>2 < h, (t>3-<t>i < 27r. (3.5) 

This is a 2-dimensional simplex, which can be described in a more symmetrical 
way as follows. Let 

ei = (1, 0), 62 = V3), 63 = - V3) (3.6) 

be the corners of an equilateral triangle in the plane. To any vector (^i, ^2? V^'s) 
in the fundamental domain, we may associate a point v in the plane through 

V = I {(piBi + (j>2e2 + (pae^) . (3.7) 

This relationship is one-to-one, since 

(j)„ = vea (3.8) 

for all a = 1,2, 3. The image of the mapping is the interior of an equilateral 
triangle with corners 

a = 0, b = — ^61, c = — 63 (3.9) 

(see fig. 1). In other words, any point v in this area determines a vector 
(^1, (f)2,(t>'i) in the fundamental domain and vice versa. 

The triangle shown in fig. 1 is mapped onto itself under refiections at the 
dashed lines. A refiection at the line passing through the point c corresponds 
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Fig. 1. Fundamental domain in the v-plane (area enclosed by the 
triangle abc). The point A corresponds to the angles (0i, 02,03) chosen 
for the definition of the running coupling. 



to a transformation 0„ where 

4>i = -<t>i - 

h = -h + y, (3.10) 

- ^ 2x 
93 = -92 + -y- 

If we define the boundary field W through eqs.(3.1),(3.2), the replacement of 
(j)a by (j)a amounts to a charge conjugation followed by a central conjugation 
(a gauge transformation which is periodic up to a central phase factor [12]). 
The transformation (3.10) thus corresponds to a discrete symmetry of the 
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Schrodinger functional and so do the other reflections at the dashed lines in 
fig. 1. 



3.3 Definition of the running coupling 

We now choose a 1-parameter curve of vectors (^i, (j)2, (ps) and (f)^, (p'^) in 
the fundamental domain. The efl"ective action 

r = -lnZ (3.11) 

then becomes a function of the curve parameter rj and a renormalized coupling 
may be deflned by difl"erentiating T with respect to rj. 
The curve of angles (^i, V^'s) that we decided to take is given by 

± ^ 1 

(pi = r)Wi - -, wi - 1, 

1 

(f>2=T]W2, W2 = --+V, (3.12) 



TT 1 



V. 



3' 2 

This is a straight line which passes through the point A at 7/ = and so is 
contained in the fundamental domain for small values of 7/ (see flg. 1). The 
direction {wi,W2,W'i) of the line may be adjusted by changing the parameter 

V. 

The other angles (f)'^ are chosen to be related to (f)^ through the symmetry 
(3.10), i.e. we set (f)'^ = (j)^. The partition function is then invariant under a 
discrete symmetry transformation of the gauge fleld, which is a combination of 
a time reflection, a charge conjugation and a central conjugation. In particular, 
the boundaries at a;" = and = L are treated equally. 

We can now deflne a renormalized coupling g^ through 



k_ 

i2' 



(3.13) 



where 

= 12(i/a)2[sin(6l) + sin(26l)], 6* = ix(a/i)2. (3.14) 

This choice of proportionality constant guarantees that g^ will be equal to g^ 
to lowest order of perturbation theory, for all values of the lattice spacing. The 
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reason for setting z/ = is that in numerical simulations of the Schrodinger 
functional the statistical errors on the coupling turn out to be minimal in this 
case. Note that the box size L is the only external scale appearing in the 
definition (3.13), i.e. is a running coupling defined at distance L. 
For general i/ we have 



dT 



= k\^^-vv\, (3.15) 



where v is another renormalized quantity. A moment of thought reveals that 
V is independent of v. Our interest in this quantity arises from the fact that 
it can be computed with little extra work and that it may be used to test the 
universality of the Schrodinger functional. 

3.4 Relation to other schemes 

In the continuum limit, and at sufficiently high energies (small L), diff"erent 
running couplings may be related in perturbation theory. The coupling defined 
above can be computed to 1-loop order using the techniques described in ref.[2]. 
Since there is no new element involved here, we only quote the results of our 
calculations. 

Let us define 

«(?) = ^, g = iA. (3.16) 

The connection between the MS scheme of dimensional regularization [13] and 
our finite volume scheme then is 

"ms = " + ^1"^ + • • •' fci = 1.25563(4), (3.17) 

where both couplings are evaluated at the same momentum q. 
Another coupling in infinite volume, 

«qq(g) = Fi^(^), q = lA, (3.18) 

is obtained from the force F(r) between static quarks at distance r. Combining 
the expansion above with the 1-loop results of refs. [14,15], one finds 

aqq = a + Ziia^ + /ii = 1.33776(4), (3.19) 

i.e. to this order there is practically no diff"erence between aqq and cij^. 
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We mention in passing that v vanishes to lowest order of perturbation 
theory. The leading term in the continuum limit is 

V = 0.0694603(1), (3.20) 

while a determination of the next order correction (a term proportional to a) 
would require a 2-loop calculation. 

3.5 Renormalization group 

We again assume that the continuum limit has been taken and define the 
Callan-Symanzik /3-function through 

m = -L§^- (3.21) 

From eq.(3.17) and the known perturbation expansion of the /3-function in 
the MS scheme, we infer that 

oo 
n=0 

with [16-19] 

6o = ll(4x)"% 6i = 102(4x)"\ (3.23) 

The 3-loop coefficient 62 depends on our choice of running coupling and is 
presently not known. 

In the following a key role is played by the step scaling function a[s,u). 
For any given scale factor s and initial value u = g^[L), the coupling u' = 
g^[sL) may be computed by integrating the renormalization group equation 
(3.21) (assuming the /3-function is known), u' is a well-determined function 
of s and u, and so we may define 

a{s,u) = u'. (3.24) 

In other words, the step scaling function is an integrated form of the (3- 
function, which tells us what happens to the coupling if the box size is changed 
by a factor s. 
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Table 1. 1-loop coefficients 6i{a/L) and ei{a/L) [eqs.(3.26),(3.28)] 



L/a 


^1 


ei 


6 


-0.00394 


0.0677 


8 


-0.00194 


0.0336 


10 


-0.00117 


0.0204 


12 


-0.00079 


0.0138 


14 


-0.00057 


0.0100 


16 


-0.00043 


0.0076 



3.6 Cutoff effects 

We finally address the question of how strongly the evolution of the coupling 
is affected by the underlying space-time lattice. Let us consider two lattices 
with size L and sL at the same bare coupling go (we assume L and sL are 
integer multiples of a). Suppose u and u' are the values of on these lattices 
and define the lattice step scaling function through 

S(5,u,a/i) = u'. (3.25) 

Close to the continuum limit one expects that the relative deviation 

Ti(2,u, a/ L) — a(2,u) ^ , ^ , o /„ „„x 

,J . ' = 6^{alL)u + 82{alL)u^ + ... (3.26) 

cryl, u) 

converges to zero with a rate roughly proportional to aj L. Since we are using 
an improved action, the 1-loop coefficient Si[a/L) is actually decreasing more 
rapidly (see table 1). The next term in the expansion (3.26) is not improved, 
however, and so, barring accidental cancellations, will be of order a/L. 

The largest coupling at which the step scaling function will be computed 
is u = 2.77. Table 1 thus suggests that it is only weakly affected by the lattice 
cutoff. Whether this is really true, must of course be verified by performing 
numerical simulations at different values of the lattice spacing (cf. subsect. 5.1). 
One should also not conclude that the lattice effects are generally small as a 
result of our use of an improved action. A sobering example here is the quantity 
V defined in subsect. 3.3. On the lattice we have 

v{L) = n{f{L),a/L), (3.27) 
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where 0(u, aj L) converges to a universal function u){u) in the continuum limit. 
As may be seen from table 1, the relative deviation 

^ ' '\ ^ = €1 a/i +^2 a/i 3.28 

is larger in this case. 



4. Numerical simulation 



The running coupling is inversely proportional to the expectation value 

of 

A5 = 1^ . (4.1) 

This observable may be computed with little effort, for any given field con- 
figuration, and so it is possible to evaluate its expectation value by numerical 
simulation of the Schrodinger functional. We here describe some of the tech- 
nical details of this calculation and list the data obtained in this way. 

^.1 Monte Carlo algorithm 

To generate a representative ensemble of gauge field configurations, simulating 
the partition function (2.2) with the required boundary conditions, we use a 
"hybrid over-relaxation" (HOR) algorithm (for a recent review and references 
see ref.[20]). The basic cycle consists of 1 heatbath update sweep through 
the lattice followed by N over-relaxation sweeps. All links are visited and 
updated once in every sweep. The program processes one time-slice after the 
other, and each time-slice is further divided into sublattices which can be 
updated independently of each other. 

Both the heatbath and the microcanonical algorithm proceed through 
embedded SU(2) subgroups [21]. We use the three obvious subgroups which 
leave one of the basis vectors in the fundamental representation of SU(3) 
invariant. The situation then is essentially the same as in the SU(2) theory 
discussed in ref.[3]. In particular, we again employ the heatbath algorithm of 
Haan and Fabricius [22,23], and the exactly energy preserving moves of the 
embedded SU(2) links are carried out in the usual way. This is technically 
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easier than proposing approximately microcanonical steps in the full SU(3) 
group followed by a non-trivial acceptance decision [24,25]. 

4.2 Autocorrelations and timing 

The number N of over-relaxation sweeps per heatbath sweep is a free param- 
eter of the HOR algorithm. In spin models and free field theory it is known 
that minimal autocorrelation times result, if iV is taken to be an approximately 
constant multiple of the correlation length in lattice units. In our computa- 
tions the relevant scale is L and N ~ L/2a proves to be a nearly optimal 
choice. We also found it profitable to determine AS" after each over-relaxation 
sweep. Little additional CPU time is required for this, while a significant loss 
of information is avoided when N is large. 

On the larger lattices most of the time is spent doing over-relaxation 
sweeps. At L/a = 20 our program achieves an average link update time of 15 
/xsec. This is about a factor 8 slower than the analogous SU(2) program, but 
since the variance of AS" is smaller here, we only need twice as much computer 
time to obtain the running coupling to the same relative accuracy. 

In most cases the observable AS" is fiuctuating about its mean with a 
nearly Gaussian distribution. Rare excursions to very small or negative val- 
ues were however observed at the largest values of the renormalized coupling 
considered. These correspond to some long-time autocorrelations and make a 
reliable error estimation difficult. Inspired by the multi- canonical technique 
for first order phase transitions [26] the problem could be overcome by a sim- 
ple modification of the sampling procedure. Further details are reported in 
appendix A. 

4.3 Simulation results 

To study the evolution of the running coupling we have simulated pairs of 
lattices with sizes L and sL at the same bare coupling. The results of the 
computations are listed in table 2. There are 5 blocks of data, corresponding 
to 5 fixed values of g^{L). In each case the bare coupling was tuned, using a 
reweighting technique, to obtain the desired value of g^{L). The scale factor 
s was set equal to 2 except for the last evolution step where s = 3/2. 

The statistical errors were calculated by jacknife binning and checked by 
summing the autocorrelation function over an appropriate time interval when 
no reweighting was necessary. The precision in g^ given in table 2 was achieved 
by accumulating around 160k sweeps on the lattices with L/a = 16 and 200k 
sweeps if L/a = 20. On a single CRAY YMP processor this corresponds to 
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Table 2. Evolution of the renormalized coupling at fixed (3 = &/g1 
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L/a 




f{sL) 


v{L) 


v{sL) 


8.7522 


5 


1.2430(12) 


1.4284(53) 


0.073(2) 


0.063(6) 


8.8997 


6 


1.2430(13) 


1.4270(50) 


0.067(2) 


0.053(6) 


9.0350 


7 


1.2430(15) 


1.4230(50) 


0.064(3) 


0.063(5) 


9.1544 


8 


1.2430(14) 


1.4250(58) 


0.071(2) 


0.065(6) 


8.1555 


5 


1 4300f211 


1 6986f711 


072f3l 


054f5l 


8.3124 


6 


1.4300(20) 


1.6859(73) 


0.069(3) 


0.064(5) 


8.4442 


7 


1.4300(18) 


1.7047(77) 


0.065(2) 


0.064(6) 


8.5598 


8 


1.4300(21) 


1.6966(90) 


0.066(2) 


0.071(8) 


( .ODO 1 









U.U ( U(^z j 


U.Udd(^o j 


7.7170 


6 


1 6950f26l 


2 091flll 


066f2l 


058f5l 


7.8521 


7 


1.6950(28) 


2.112(10) 


0.064(2) 


0.052(5) 


7.9741 


8 


1.6950(28) 


2.096(11) 


0.063(2) 


0.060(5) 


8.1650 


10 


1.6950(31) 


2.092(15) 


0.066(2) 


0.053(6) 


6.9671 


5 


2.1000(41) 


2.750(17) 


0.066(2) 


0.050(5) 


7.1214 


6 


2.1000(39) 


2.783(18) 


0.059(2) 


0.047(4) 


7.2549 


7 


2.1000(43) 


2.774(16) 


0.058(2) 


0.059(4) 


7.3632 


8 


2.1000(45) 


2.772(17) 


0.059(2) 


0.059(4) 


7.5525 


10 


2.1000(42) 


2.743(24) 


0.058(2) 


0.040(6) 


6.5512 


6 


2.7700(69) 


3.489(22) 


0.053(2) 


0.047(3) 


6.7860 


8 


2.7700(73) 


3.448(28) 


0.051(2) 


0.043(3) 


6.9748 


10 


2.7700(76) 


3.487(29) 


0.052(2) 


0.045(4) 


7.1190 


12 


2.770(11) 


3.496(42) 


0.051(3) 


0.050(5) 


s = 1 at g2 


(i) = 


2.770 and s = 2 in 


all other cases 







180 and 600 hours of CPU time respectively. 

An additional set of simulations were performed to determine the bare 
coupling as a function of L/a at fixed g^{L) = 3.48 (table 3). These results 
will be used in subsect. 5.3 to make contact with the physical low-energy scales 
of the theory. 
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Table 3. Bare couplings vs. lattice size at f{L) = 3.48 



L/a 


/3 


L/a 


/3 


4 


5.9044(38) 


8 


6.4527(46) 


5 


6.0829(35) 


9 


6.5539(80) 


6 


6.2204(25) 


12 


6.775(13) 


7 


6.3443(43) 


15 


6.973(10) 



5. Discussion of results 



5.1 Extrapolation to the continuum limit 

From the numbers listed in table 2 one obtains the lattice step scaling function 

^{s,g\L),a/L) = f{sL) (5.1) 

at 5 values of g^{L) and for a range of a/L. We now pass to the continuum 
limit a/i = by extrapolating these data, using an ansatz of the form 

S(s, -g\ a/L) = a{s, f) {l + p{s, f) a/L] . (5.2) 

As shown by fig. 2 the fit works perfectly. Within errors there is no observable 
cutoff" dependence, i.e. the slope p[s,g^) is compatible with zero in all cases. 
The results of the fit will be discussed below. 

Compared to the SU(2) theory the lattice step scaling function here is 
much less dependent on the cutoff. Presumably this is due to our use of an 
improved action. It is interesting to note in this connection that the SU(2) 
data of ref.[3] can be fitted over the whole range of couplings with an error 
term of the form p{2,g^) = p2g^ ■ Moreover, the fit parameter p2 is found to 
have the same sign and order of magnitude as the coefficient of the a/L lattice 
correction which one obtains at 1-loop order of perturbation theory (recall that 
in the SU(2) theory the simulations were performed with the Wilson action). 
So it is quite plausible that the leading cutoff effects may be cancelled by a 
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3.5 



2.5 



1.5 



S(3/2,2.770,a/i) 
^ . ^ 



S(2, 2.100, a/i) 



Y 




m 




S(2, 1.695, a/i) 
• • • • 


• — 












S(2, 1.430, a/i) 


• — 


G 


• • • 

• • • 


• — 



J I I L 



S(2, 1.243, a/i) 

J I I \ I I I I I 



J I L 



J I L 



0.05 0.1 0.15 

a/L 



0.2 



0.25 



Fig. 2. Extrapolation of the lattice step scaling function to the con- 
tinuum limit. The points at a/L = represent the extrapolated values. 

perturbative improvement of the action. Further studies are however needed 
to completely clarify the issue. 

We now proceed to discuss the cutoff dependence of v, or, more precisely, 
of the function a/L) [eq.(3.27)]. The data listed in the last two columns 

of table 2 show that v is hardly changing as a/L is made smaller at fixed g^. 
Note that in most cases we have results for a/L from 0.2 down to 0.05 (small 
variations of can here be ignored, because is a fiat function of g^). This 
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Table 4. Values of the step scaling function 





s 






1.243 


2 


1.416(16) 


1.428 


1.430 


2 


1.703(24) 


1.684 


1.695 


2 


2.095(25) 


2.071 


2.100 


2 


2.771(41) 


2.732 


2.770 


3/2 


3.474(65) 


3.397 



nicely confirms the universality of the Schrodinger functional [2], although it 
must be said that the test is not too significant, since the statistical errors are 
quite large. 

For the extrapolation to the continuum limit, we use the formula 

n{g\ alL) = uj{f) {1 + ei(a/i) + a/i} , (5.3) 

where ei[a/L) is the exact 1-loop lattice correction (cf. subsect. 3.6). The 
remaining cutoff" effects (those proportional to A) then are compatible with 
zero and certainly less than 10% in the whole fit range. 

5.2 Evolution of the running coupling 

Our results for the step scaling function in the continuum limit are compiled 
in table 4. Up to small mismatches the couplings are arranged so that one can 
step down the energy scale successively, starting from the smallest coupling, 
g^ = 1.243, and taking a zigzag course through the table until one arrives at 
the largest coupling. The total scale factor then is approximately equal to 24. 

Following ref.[3] the integration of the step scaling function can be per- 
formed with all errors and mismatches properly accounted for. At this point 
it is convenient to use imax, the box size at which f = 3.48, as a reference 
scale. The results are listed in table 5. 

It is now interesting to compare the evolution of the coupling with the 
predictions of perturbation theory. In the last column of table 4 we have in- 
cluded the values of the step scaling function which one obtains by integrating 
the renormalization group equation (3.21) using the 2-loop approximation to 
the /3-function. There is hardly any difference between these numbers and 
the simulation results. Only at the larger couplings does one see a small sys- 
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Table 5. Running coupling at scales L given in units of L 



L 1 Lpaax 








1.000 


3.480 


0.165(9) 


1.695 


0.664(19) 


2.770 


0.084(6) 


1.430 


0.332(14) 


2.100 


0.040(4) 


1.243 



tematic deviation which is statistically significant, because the errors on the 
step scaling function are uncorrelated (they are extracted from disjoint sets of 
simulations). 

In view of these remarks it is not surprising that a perfect representation 
of the data listed in table 5 can be obtained by integrating the renormalization 
group equation (3.21), starting from = 3.48 and assuming 

m = -bof - hrf - bff. (5.4) 

A correlated least squares fit yields 

bf = 1.5(8) X (4x)-3 (5.5) 

with a per degree of freedom equal to 1.4/4. There is no reason to expect 
that coincides with 62. But it is interesting to note that is about twice 
as large as the exact 3-loop coefficient in the minimal subtraction scheme of 
dimensional regularization [27]. In particular, the correction appears to be 
reasonably small, when the /3-function is written as power series in a. 

In the range of couplings considered, there is, therefore, no discernible dis- 
crepancy between perturbation theory and the simulation results. We mention 
in this connection that the function ijo{g^) has a perturbation expansion of the 
form 

u;{f) = u;o+u;ig'^ + ..., (5.6) 

where the 1-loop coefficient uiq is given by eq.(3.20). u)i is presently not known, 
but the simulation data can be fitted rather well with eq.(5.6) and a reasonable 
eff"ective 2-loop coefficient (see fig. 3). 

We would like to emphasize, however, that our results do not prove that 
perturbation theory provides a good approximation to all quantities of interest 
up to couplings as large as 3.48. Such a general statement is bound to be false 
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Fig. 3. Plot of V = w[g^). The dashed line is a fit, using eq.(5.6) with 
an adjustable 2-loop coefficient, to all data points except the one at the 
largest coupling. 

and the running coupling in our scheme may very well turn out to be an 
exceptional case. 

5.3 Physical units 

We now need to relate imax to the dynamical low-energy scales of the theory 
in infinite volume. As discussed in sect. 1, we decided to take the distance t-q 
as the basic reference scale. Our first goal is thus to compute the conversion 

factor imax/f'O- 

The heavy quark force F[r) has previously been calculated through nu- 
merical simulation of the Wilson theory on large lattices. Here we use the 
results of refs.[8,28] to evaluate r^/a for various values of the bare coupling 
(see table 6). The details of the computation are as in the SU(2) theory [5]. 
In particular, the errors quoted include an estimate of the uncertainty which 
arises from a possible admixture of excited energy levels to the heavy quark 
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Table 6. Low-energy scales at various values of /3 



13 


ro/a 


-^max/^ 


Lmax/ fO 




roVK 


reference 


6.0 


5.44(26) 


4.516(13) 


0.830(40) 


0.0513(25) 


1.23(7) 


[10,28] 


6.2 


7.38(25) 


5.815(11) 


0.788(27) 


0.0250(4) 


1.16(3) 


[10,28] 


6.4 


9.90(54) 


7.488(21) 


0.756(42) 


0.0139(4) 


1.16(6) 


[10,28] 


6.5 


11.23(21) 


8.498(31) 


0.756(14) 


0.0114(4) 


1.20(3) 


[8] 



F{r) and the value of K are taken from the references quoted in the last column 



potential. 

The values of L^g^^/a listed in the third column of table 6 have been 
obtained from table 3 using the interpolation formula (3 = ai +a2 ln(i/a). It is 
now trivial to form the ratio L^g^^/ro (4th column of table 6). L^i^^/ro depends 
on the cutoff and must be extrapolated to the continuum limit. Assuming a 
linear dependence on the lattice spacing, the extrapolation yields 

Lm../ro = 0.674(50) (5.7) 

at a/L^i^x = (see fig. 5). The data are certainly consistent with this proce- 
dure, but the errors are quite large and there are only 4 data points. To be 
able to do better we would need a more precise determination of the heavy 
quark force in the range of /3's considered. One more value of r^/a at say 
l3 = 6.8 would further stabilize the extrapolation (the data of ref.[10,28] are 
from a physically small lattice and do not reach sufficiently large distances r). 
Such calculations are possible with the available hardware and will hopefully 
be done at some point. 

In table 6 we have included the values of the string tension quoted in 
refs.[8,10]. As an additional check on the systematic uncertainty arising from 
the extrapolation of F[r) to large distances, we have also made our own fits. 
In particular, the fit range was varied and an r""^ correction besides the usual 

term was included in the fit formula. Taking into account the spread 
of values obtained in this way, we get rQ\/K = 1.22(8), the result quoted in 
sect. 1. 

For the purpose of illustration and to further the intuitive understanding 
of our results, we now convert to more physical units by setting t-q = 0.5 fm. 
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Fig. 4. Extrapolation of Lma^/ro to the continuum limit. The left- 
most point represents the extrapolated value (5.7). 

We then deduce that imax = 0.337(25) fm, and the smallest box size cov- 
ered by table 5 is hence roughly equal to 0.014 fm. For the string tension we 
get = 481(32) Me V, which is somewhat larger than the values usually 

quoted. These are obtained from the phenomenological heavy quark potential 
at distances around 1 fm and also from other sources, such as the relativistic 
string model and the observed slope of the p meson Regge trajectory. In the 
case of the heavy quark potential, the difference arises from the fact that the 
theoretical potential (as computed in lattice gauge theory) is steeper than the 
phenomenological potentials commonly used to fit the charmonium spectrum 
[5]. 

In fig. 5 the running coupling a[q) [eq.(3.16)] is plotted as a function of 
the momentum q given in physical units. The error bars on the data points in 
this figure are barely visible. They represent the statistical errors as inferred 
from table 5, but not the overall scale uncertainty which originates from the 
conversion factor (5.7). The latter amounts to a multiplication of the energy 
scale by a constant factor and so has no bearing on the scaling properties of 
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1 \ — \ — \ — [" 




1 10 

q[GeV] 

Fig. 5. Comparison of numerically computed values of the running 
coupling (data points) with perturbation theory. The full line represents 
the fit described in subsect. 5.2, while the dashed (dotted) curve is obtained 
by integrating eq.(3.21), starting at the right-most point and using the 2- 
loop (1-loop) /3-function. 

the coupling. 
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5-4 Computation of 

Using the fit of our data discussed in subsect. 5.2, and the conversion factor 
(5.7), the running coupling ot(q) is obtained at all momenta q between 1.5 X r^^ 
and 37 X r^^ . In view of fig. 5 we might actually go to even larger momenta 
with little risk of running into extrapolation errors. 

We may now convert to the MS scheme of dimensional regularization 
by applying the 1-loop formula (3.17). Of course this is only possible at 
sufficiently high energies where the coupling is small. A typical result is 

aj^(g) = 0.1108(23)(10) at g = 37 X t-q-^ (5.8) 

where the first error is the combined statistical error from the fit and the 
conversion factor (5.7), while the second is equal to a[q)'^ and thus indicates 
the expected order of magnitude of the 2-loop correction in the matching 
formula (3.17). This latter error cannot be reliably assessed until the necessary 
2-loop computations have been performed. But it is reassuring to note that a 
correction of the quoted order of magnitude will result if 62 is in the range of 
the effective 3-loop coefficient [eq.(5.5)]. 

5.5 Comparison with other computations of the running coupling 

In the work of El-Khadra et al. [6] the running coupling is obtained in two 
steps. One first determines the mass splitting Am between the IP and IS 
(quenched) charmonium levels on a large lattice, where finite volume effects 
can be neglected. If one takes Am as the low-energy reference scale, this 
computation may be regarded as giving the lattice spacing a in physical units. 
The running coupling at the momentum q = x/a, which is now known in 
physical units, is then calculated using the "improved" 1-loop formula [29-32] 

a^^(g) = ao/(l + 0.309 X ao). (5.9) 

The modified bare coupling cto occurring here is defined through [29] 

ao = ffoV(47rP), P = i(trC/(p)), (5.10) 

and the plaquette expectation value P is to be computed at the given value 
of fl'o- 

To be able to compare with our results in a more direct way, we use r^ 
instead of Am as the reference scale. From the second column of table 6 
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we can then read off the lattice spacing in physical units for a range of bare 
couplings. So if we choose (3 = 6.5, for example, we have q = 35.3(7) X r^^ 
and, noting P = 0.6384 [8], one obtains a-g^i^q) = 0.1111. Our value for the 
MS coupling at this momentum is 0.1119(23)(10) and a similar agreement is 
found at the lower values of (3 in table 6. 

The corrections to eq.(5.9) are thus at most a few percent in this range 
of bare couplings. Note that if one employs another reference scale, such as 
Am, the statement remains true provided the corresponding conversion factor 
can be shown to be independent of the lattice spacing, to a sufficient degree 
of precison, at the values of (3 considered. 

In refs.[7-ll] the coupling aqq(g) [eq.(3.18)] is computed from lattice data 
on the heavy quark force. A subtle fit and subtraction of cutoff effects is needed 
at short distances, where the raw data are strongly dependent on the lattice 
action employed. The values of aqq obtained in this way are significantly larger 
than those one gets by combining the 1-loop formula (3.19) with our results 
for a. At /3 = 6.5 and r ja = 2.5322, for example, the value quoted in ref.[8] is 
0.248(2)(1), while our result is 0.205(7)(5). There is no reason to be worried 
about this discrepancy, because the distances at which the comparison can 
be made are rather large in physical units. So it could well be that the non- 
perturbative corrections to the matching formula (3.19) are not negligible. It 
is also not evident to us that the cutoff effects are indeed under control, when 
ttqq is determined from the heavy quark potential at distances as low as 2 or 
3 lattice spacings. 



6. Concluding remarks 

Compared to our earlier study of the SU(2) Yang-Mills theory [3], the re- 
sults obtained here are very similar. In particular, we again find that the per- 
turbative scaling regime and the low-energy domain of the theory are smoothly 
connected, with no complicated transition region. We have also been able to 
determine a-g^i^q) at momenta q up to about 14 GeV, with an estimated pre- 
cison of a few percent. An interesting observation is that the relation 

2«Ms(«)|sU(2) = 3"ms(«)|sU(3) (6-1) 

holds within errors, if q is given in units of tq. Eq.(6.1) is consistent with 
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the evolution of the coupling up to 3-loop order of perturbation theory and 
suggests that Na-^ is only weakly dependent on the number N of colours. 

A critical step in our computations is the evaluation of the conversion 
factor L^i^x/ro, which is needed to make contact with the low-energy physics 
of the theory in infinite volume. The situation is more difficult in the SU(3) 
theory than in the SU(2) case, because the large volume lattices which one has 
been able to simulate to date have significantly larger lattice spacings. More- 
over, the statistical errors on the heavy quark force at intermediate distances 
(where t-q is determined) are far from negligible. We hope that these deficits 
can be removed in the near future. 

The finite-size scaling technique employed in this paper is expected to be 
useful in QCD, too. Besides the running coupling one is here also interested 
to determine the renormalized quark masses. The technical difficulties are 
considerable, however, and it may take a number of years until a plot like fig. 5 
can be drawn for QCD with a non-zero number of light quarks. In particular, 
the renormalizability of the Schrodinger functional must be discussed and the 
performance of numerical simulation algorithms for dynamical fermions on 
physically small lattices must be evaluated. 

We are indebted to A. Kronfeld for correspondence and to G. Bali for 
communicating simulation results on the heavy quark force prior to publica- 
tion. The computations have been performed on the CRAY computers at 
HLRZ and CERN. We thank these institutions for their support. 



Appendix A 

At the largest couplings considered [g^ ~ 3.48), a straightforward simula- 
tion of the Schrodinger functional is unsatisfactory, because some rare fiuctua- 
tions of the observable AS" lead to uncomfortably large autocorrelation times. 
We here describe how one can get around this problem by using a variant of 
the multi- canonical technique. 

The history and distribution of AS" for a "normal" simulation run are 
shown in figs. 6a and 6b. The corresponding data for a run at f = 3.45 
are displayed in figs. 6c and 6d (cf. table 2). In both cases we plotted block 
averages from 81 successive sweeps as individual "events" to smooth out the 
short time noise. To a good approximation they can be treated as statistically 
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Fig. 6. Histories and distributions of AS at L/a = 12. a,b: /3 = 
7.1214; c,d: /3 = 6.7860; e: the same, with modified sampling (7 = 0.05), 
and f: after reweighting. 

independent in normal runs. Spikes as the one clearly seen in fig. 6c however 
lead to longer autocorrelation times. They are associated with a long tail 
of the distribution of AS" at low values of AS" (cf. fig. 6d). The problem is 
that the spikes influence the expectation value of AS" significantly, but are not 
sufficiently frequent for a reliable determination of the statistical error. 

The statistics in the tail of the distribution of AS" can be enhanced by 
sampling the gauge fields with the modified Boltzmann weight 



P[U] oc exp{-S[U] - jASiU]). 



(A.l) 



The extra term proportional to 7 is easily accounted for in the simulation 
algorithm. For 7 > the simulation is now biased towards smaller values of 
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AS" and the desired expectation value is obtained through 



(AS) = —. r!'"^ . (A.2) 

(exp(7A5))~ 

The optimal value of 7 depends on the observed width of the distribution 
of AS" and requires some tuning. For 7 = 0.05 the distributions of AS" and 
AS'exp(7AS')/(exp(7AS'))~ in the P-ensemble are as shown in figs. 6e and 
6f. Evidently the problem has been overcome in this way. 
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